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(2J[)| We look in detail at the process of mapping an astrophysical initial model from 

a stellar evolution code onto the computational grid of an explicit, Godunov type 
code while maintaining hydrostatic equilibrium. This mapping process is com- 
mon in astrophysical simulations, when it is necessary to follow short-timescale 
dynamics after a period of long timescale buildup. We look at the effects of spatial 
resolution, boundary conditions, the treatment of the gravitational source terms 
in the hydrodynamics solver, and the initialization process itself. We conclude 
with a summary detailing the mapping process that yields the lowest ambient 
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^ '_ 1. Introduction 



velocities in the mapped model. 

Subject headings: hydrodynamics, methods: numerical, stellar dynamics 



Many astrophysical phenomena involve a dramatic change between timescales of interest 
- the slow convection and simmering in the interior of a white dwarf followed by ignition 
of a Type la supernova, for example, or the accretion of a layer of fuel onto a white dwarf 
or neutron star leading to ignition and runaway at the base of the layer, producing a nova 
or an X-ray burst. These two regimes are difficult to follow with a single hydrodynamic 
algorithm because of the disparity of the relevant timescales. Modeling these long timescale 
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events requires an implicit or anelastic hydro dynamic method; short timescale events require 
explicit hydrodynamic methods that can capture the transient phenomena. Often, a one- 
dimensional stellar hydrodynamics code is used to follow the accretion process until just 
before the fuel reaches the ignition temperature. Simulations of the flash following the 
ignition requires a multidimensional hydrodynamics code. Matching the two regimes is a 
difficult process, and can introduce numerous errors into a calculation. 

Simulations in the atmosphere or interior of a star or compact object often begin with 
an initial model that was generated by a 1-d implicit code (see for example Weaver, Zim- 
merman, & Woosley 1978; Glasner, Livne, & Truran 1997; Sugimoto et al. 1982; Young et 
al. 2001; Starrfield et al. 1998), that evolves the system through the long timescale processes 
(accretion, slow convection, simmering of nuclear fuel) until just before the short timescale 
dynamics begin. This 1-d initial model is then used as input to a multidimensional, explicit 
hydrodynamics code (see for example Glasner, Livne, & Truran 1997; Kercek, Hillebrandt, 
Truran 1998; Zingale et al. 2001; Bazan & Arnett 1998; Kane et al. 2000). The mapping 
of a 1-d hydrostatic initial model onto a multidimensional grid is the focus of the present 
paper. In the absence of any perturbations or external forces, the system should remain in 
hydrostatic equilibrium after the mapping. 

The mapping process can introduce a variety of errors. It is common for the two codes 
to use different equations of state (EOS), which can have a large effect on the structure 
of the atmospheres. Even if the basic physical components are the same, the treatment 
of physical details (for example, Coulomb corrections or ionization) may differ, leading to 
differences in the pressure of a fluid element between the two codes, even with the same 
density, temperature, and composition. Simply updating the thermodynamics of the initial 
model with the new EOS will most likely push it out of hydrostatic equilibrium. One- 
dimensional initial models are almost always created using mixing length theory to describe 
the convection, implicitly assuming a velocity field necessary to transport the energy as 
required. If the 1-d model was convectively unstable, it is unclear how to define the 2- or 
3-d velocity field consistent with the 1-d input velocity field. When mapping into higher 
dimensions, there is not enough information to set the velocities properly. Typically the 
velocities are zeroed during this transition. 

Finally, it is unusual for the number of points and the grid spacing in the initial model to 
match that in the multidimensional grid. The initial model may have come from a Lagrangian 
code and will need to be converted into an Eulerian coordinate system via 



The initial model will then need to be interpolated onto the new grid, which will introduce 
even more hydrostatic equilibrium errors. 



dm = 4iir 2 p(r)dr. 
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Once the model is on the new grid, differences in the hydro dynamical algorithms can 
cause problems. Although the two codes are solving the same equations, it is also very likely 
that the discretization used in the two codes is different. The very definition of the quantities 
on the grid may also differ; the mapping may proceed from a 1-d implicit finite-difference 
code, where the values on the grid may represent nodal points, to a multidimensional finite- 
volume code, where the values represent cell averages. If everything else were constant 
between the codes, the differences in the discretization and the definition of the variables 
(e.g. pointwise values vs. cell averages) is enough to upset the hydrostatic equilibrium. 
Poor boundary conditions can drive velocities on the grid, pushing our initial model out 
of hydrostatic equilibrium. Ideally, the boundaries should match the physics of the initial 
model and present a smooth state to the hydrodynamics solver. 

In this article, we look at bringing an initial model to hydrostatic equilibrium by study- 
ing the effects of the initialization method, the boundary conditions, and the solver itself. An 
initial model atmosphere is mapped onto the computational grid, transverse to the direction 
of gravity. Since we are not concerned with perturbing the model after the mapping, all the 
calculations presented here will be one-dimensional for computational efficiency. We con- 
sider three different initial models — a simple analytic model, a simple hydrostatic equilibrium 
model with a complicated EOS, and an initial model from an implicit stellar hydrodynamics 
code. These three different models will allow us to isolate the importance of the different 
parts of the code (initialization, boundary conditions, and solver itself). Our goal is to find 
an optimal configuration that allows us to hold an initial model from a different hydrody- 
namics code in equilibrium in the present code until a perturbation or external force that 
we impose disturbs it. Furthermore, we want the resulting mapped hydrostatic model to 
be as close to the original initial model as possible. Regions where disturbances have not 
yet propagated should remain quiescent for tens or hundreds of dynamical timescales. The 
absolute magnitude of the velocity that a simulation can tolerate will be problem dependent; 
however, it must be considerably less than the velocity of any dynamics we wish to study 
(e.g., burning front speed, convective speed, and certainly the sound speed.) 

For the hydrodynamics, we chose to use PPM (Colella & Woodward 1984), a widely used 
Godunov method (Godunov 1959). PPM solves the Euler equations in conservative form, 
using a finite-volume discretization, guaranteeing conservation. PPM is a shock-capturing 
scheme, which is desirable for the modeling of the rapid transients in the short-timescale 
regimes we wish to study. The implementation of PPM we use is contained in the FLASH 
code (Fryxell et al. 2000), and is based on PROMETHEUS (Fryxell, Muller, & Arnett 1989). 
While we use PPM to demonstrate results in this paper, only the discussion in §2.2, where we 
discuss extensions to the method to better treat the gravitational source term, is particular 
to PPM. Discussion of initialization and boundary conditions should be relevant to any 
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Godunov-type method, and perhaps other finite-volume methods. Although FLASH can 
use an adaptive mesh, in this study we run all the simulations on a uniform grid. 

This paper is organized as follows: §2 discusses the hydrodynamics algorithm employed 
in this study and improvements made to better follow a hydrostatic atmosphere. In §3 
we look at the initial models that will form the basis of our tests. §4 discusses the different 
boundary conditions considered. In §5, we show the results of a grid of calculations of each of 
the different initial models, varying spatial resolution, boundary conditions, and the details 
of the hydrodynamics. Finally, we conclude in §6. 



2. Hydrodynamics 

Understanding how a simulation code treats the hydrodynamics is critical to being able 
to accurately initialize the grid and maintain an atmosphere in hydrostatic equilibrium. The 
Godunov-type methods we consider here solve the compressible Euler equations of continuity, 

| + V.pv = 0, (2) 



momentum, 



energy, 



1 V • pvv + VP = pg , (3) 



and species advection, 



dt 
dpE 

d P X, 



+ V • (pE + P) v = pv ■ g , (4) 

, V • pA> = , (5) 

where p is the total mass density, v is the velocity, P is the pressure, E is the total specific 
energy, and Xi is the mass fraction of species i. In all discussion below, unless otherwise 
noted, it is assumed that the gravitational acceleration, g, is spatially constant and negative, 
and in the z direction so that g = gz. 



2.1. Cell-averaged quantities 

The Euler equations are expressed above in their conservation- law form, 

^ + V ■ F(U) = S(U) , (6) 
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where U = (p, pv, pE) T are the conserved quantities. Equation (6) is the differential form of 
the conservation law. Godunov methods, and other finite-volume methods, solve Eq. (6) in 
integral form: 

^/{f + v-F(U)}^ = i/ S (U) < «/, (7 > 

where V is the cell volume. Using (U) to denote the cell average of quantity U, we can 
rewrite this as 

^ + iyF(U)n.dS = (S(U)> . (8) 

The integral form has the advantage of being able to deal with discontinuities (e.g. shocks) 
because it requires less smoothness of the data (Toro 1999). The second term in the ex- 
pression above is simply the sum of the fluxes through the boundaries of the cell. Since 
whatever flux leaves one cell enters the neighboring cell, any finite-volume discretization of 
conservative equations will itself be conservative to roundoff error. 

Much of what we will need to do throughout the rest of the paper to bring a model into 
hydrostatic equilibrium will involve interpolating or extrapolating from data to the mesh. 
However, because finite-volume methods deal with cell-averaged quantities, rather than cell- 
center or nodal quantities, one cannot use usual interpolants when creating functions that 
reconstruct the discrete data. Rather than using a function f(z) such that / (zj) = fi, one 
must create a function such that (1/Sz) f*_ f(z')dz = (f) i (see Figure 1). For a uniform 
mesh in one dimension, polynomial reconstruction functions with z = taken to be the cell 
center containing the value (/) (see also Laney 1998) are for first-order: 

m = if)+1 ~ {f)o z+(f) , o) 

for second order: 

f( ^ _ (/>+! ~ 2 (/>0 + (/>-! 2 , (/>+! ~ (/>-! , - (/>+! + 26 (/>0 ~ (/>-! (m , 

f[ } ~ 25* Z + 25 Z+ 24 ' 1 j 



for third order: 



f{z) = Az 3 + Bz 2 + Cz + D , (11) 
(/> + i - 3 </>„ + 3 </>_!- </>_ 2 



.4 
B 
C 
D 



65 3 

(f) +l - 2 (f) + (/)_! 
25 2 

7(/> +1 + 15(/) -27(/)_ 1 + 5(/>_ 2 
245 

- (/> +1 + 26 (/> -(/>_! 
24 
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and, finally, for a fourth order polynomial: 



/(*) = 


Az 4 + Bz 3 + Cz 2 + Dz + E , 


A = 


(/) +1 -4(/) + 6(/)_ 1 -4(/)_ 2 + (/)_3 


245 4 


B = 


3 </) +1 - 10 (f) + 12 (/)_, - 6 (f)_ 2 + (/)_ 3 


125 3 


C = 


7(/>+i-12</>o + 2</)_ 1 + 4(/)_ 2 -</)_3 


165 2 


D = 


9 (f) +l + 50 </>„ - 84 + 30 (/>_ 2 - 5 (/>_ 3 


485 


E = 


-71 + 2044 (f) Q - 26 (/)_! - 36 (/)_ 2 + 9 (/> 


1920 



These reconstruction functions will be used throughout the paper. 

A high-order Godunov method solves the Euler equations in conservation-law form using 
Eq. (8). In one approach, the first step in the calculation of fluxes between the finite- volume 
zones is to perform a reconstruction to represent values of the variables continuously through 
the zone. In PPM, this reconstruction step is particularly involved. Care is taken to introduce 
no new maxima/minima in the polynomial representation of the solution. Furthermore, 
near shocks and discontinuities, the profile is flattened, so as to avoid oscillations. The 
reconstructed profiles are then averaged over the region that will be 'seen' by waves during 
the next timestep to generate average values on the left and right side of the interface, to 
provide left and right states as input to the Riemann problem. The Riemann solver then 
constructs the fluxes through the boundary, which are then used to update the cell-averaged 
solution vector, (U). The Riemann solver we used is based on Colella & Glaz (1985), and 
can handle an arbitrary EOS. 



2.2. Hydrostatic equilibrium 

We note that in the absence of any initial velocities, Eqs. (2) to (5) reduce to 

VP-„. ^0, ^0. , 13 ) 

These equations represent the condition of hydrostatic equilibrium (HSE). Simulations begin- 
ning from initial models satisfying VP = pg, with zero velocity everywhere should maintain 
that profile with time. If this expression is not satisfied perfectly, as differenced in the hydro- 
dynamics algorithm, an acceleration will result. The non-linear character of these equations 
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means the resulting velocities are likely to be amplified and disturb an initially hydrostatic 
model. 

In the operator-split formulation that FLASH and many other codes use, the hydro- 
dynamics and the gravity operators are not very closely coupled. In this case, maintaining 
HSE relies on the exact cancellation of two possibly large terms, VP and pg, calculated 
in two different ways. Since this cancellation won't in general be exact, spurious velocities 
result. In Godunov-type solvers, this cancellation must happen in two places: in the highly 
nonlinear Riemann solve across each cell interface during the flux calculation, and in the 
subsequent update in each zone. It is the Riemann solve which is the most problematic, and 
which we address here. 

For simulations where one of the acceleration terms dominates, the resulting small errors 
may not be significant. For nearly- hydrostatic problems, however, the spurious velocities 
could be disastrous. Thus, for these problems, one would like to forgo the splitting of 
the imposition of the gravitational and hydrodynamic acceleration, and instead include the 
gravitational effects directly in the hydrodynamic solve. Other authors, such as LeVeque 
(1998a), have suggested ways of informing Godunov-type hydro solvers of the gravity source 
term by notionally putting a constant jump into the energy and density field given to the 
Riemann solver to cancel out the source term while maintaining the cell-average quantities. 
This works quite well for both the Godunov method, where we have implemented the method 
in the FLASH code, (see Figure 2) and for a large class of higher-order Godunov methods 
(LeVeque 1998b). 

PPM, and other 'Reconstruct- Solve- Average' methods (LeVeque 1998a), use a large 
stencil to compute a smooth reconstruction to the fields on either side of the interface, and 
average over that reconstruction to generate the left and right states to the Riemann problem. 
Since the range over which the reconstruction averaged over is determined by the charac- 
teristics and isn't known a priori, it is not obvious how to nicely and efficiently compute 
a constant jump in the states to account for the gravitational source term. Nonetheless, a 
smooth analog is possible. 

In our 'modified states' version of PPM, when calculating the left and right states for 
input to the Riemann solver, we locally subtract off from the pressure field the pressure that is 
locally supporting the atmosphere against gravity; this pressure is unavailable for generating 
waves. At each cell interface (z = z i+ i/ 2 ), we subtract the pressure locally supporting the 
atmosphere by computing a reconstruction to the quantity (pg) from p(z) and g(z), and 
defining a 'wave-generating' pressure by subtracting off the integral of this quantity. In the 
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zone to the left of the interface (z = 2^+1/2), the modified pressure is 

Pftz) = Pi{z) - / (pg)i(z)dz Zi-1/2 < z < z i+1/2 . (14) 

J z 1 

In the zone to the right of the interface, the modified pressure is 

P' i+1 (z) = P i+1 (z) - (pg) i+1 (z')dz' z i+1/2 < z < z i+3/2 . (15) 



We may do this because the absolute value of the pressure doesn't matter for gener- 
ating waves, and only the pressure in excess of hydrostatic equilibrium will create motion. 
However, there are other effects where the absolute value of pressure is important — in par- 
ticular, we use P(z), not P'(z), in calculating sound speeds. We also see that, if the original 
reconstructions of pressure were continuous at the interface (Pi(z i+ i/ 2 ) = Pi+i(Zi + i/ 2 )), then 
the new reconstructions will also be. 

The reconstruction of (pg) is calculated at each point by calculating the product of p 
and g and performing the same reconstruction on this quantity as is used for quantities such 
as pressure and density. Once this (usually polynomial) reconstruction is done, integrating 
it within the cell is straightforward. The approach is shown schematically in Figures 3 and 
4, using the method described below. 

Recall that in PPM, the reconstruction is done locally by functions of the following 
form: 

/(a) = /ii + a[A/ i + / 6i (l-a)]; a= - * H , (16) 

f{P) =fni-P [A/* -P)]\ /?= *' + *_~* , (17) 

z i+\ z i-\ 

where /^j, fm, A/j, and f &i , in the notation of Fryxell et al. (2000), are the coefficients 
of the reconstruction polynomial through zone i in the normalized coordinates a or (3. The 
coefficients are chosen by the PPM algorithm to reconstruct the data consistent with other 
constraints such as monotonicity and shock capturing. The first form of the reconstruction 
is designed to be convenient when integrating through the zone from the left (e.g., when the 
zone is on the right of the interface), and the second when integrating from the right (e.g., 
when the zone is on the left of the interface). 

Given the above reconstructions, the integral of (pg) in cell to the right of the interface 
is (using Eq. 16): 

/ (pg) i+1 dz = 5z i+1 / (pg) i+1 (a)da (18) 

J z. , 1 J a=0 

! +2 
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[{pg) L i+i + a (A(pg) i+1 + (pg) 6 i+1 (1 - a))} da 



a=0 



-g(pff)e i+i« 3 + 2 ( A (P^)i+i + (PfiOe i+i) « 2 + {pg)hi+ia 



and in the cell to the left of the interface (using Eq. 17), 

rj8 



(pg)idz 



(pgWW 



0=0 



[{pg) Ri - /? (A( W ) 4 - ( W ) 6 i (1 - /?))] d/? 



0=0 

1 



1 



3 0<?)e i(3 3 + - (-A(pg)i + 0<?) 6 i) /5 2 + 



(19) 

(20) 
(21) 



(22) 

(23) 
(24) 



Note that, by construction, there are no terms in the integral constant in a or (3. 



The reconstructed pressure in zones i and i + 1 when integrating from the right and 
from the left of the interface at Zi+1/2 is 



Pi(P) = P Ri -/3[APi-P 6i (l-l3)] , 
P i+ i(a) = P Li+1 + a[AP i+1 + P 6i+1 (l-a)} . 

The new reconstructions of the pressure used in the modified states become 



pm = p Rl - p [ap; - p^- (i - p)} -/5 3 p; lbic , 



p' i+ M 



P Li+1 + a[AP! +1 + P£ +1 (l 



a 



^ -^"cubic i 



where 



AP' 



r % i 



P' 

cubic i 



AP - A Zl Upg) Li + - {A(pg)i + (pg) 6i ] 
P& i + ^Szi (A(pg)i ± (pg) 6 4 ) , 
-6zi(pg) 6i ■ 



(25) 
(26) 



(27) 
(28) 



(29) 

(30) 
(31) 



A diagram showing how this works at an interface is shown in Figure 5. 



Since these reconstructions are used to create the states for input into the Riemann 
problem, the effect of gravitational acceleration is incorporated directly into the hydrody- 
namics. Thus, its effects need not be added in later in the construction of the input states 
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to the Riemann solver; in particular, the terms corresponding to those involving gravity in 
(3^ in Eq. (3.7) of Colella & Woodward (1984) are removed. It is still necessary however to 
add the gravitational terms to the momentum and energy when doing the zone updates, as 
specified in Eq. (3.8). For example, the momentum is updated according to 

(P)T 1 « +1 = (P) n i K + H (Fi-1/2 - F i+1/2 ) + At fo)? + (p): +1 (g)^ 1 ) . (32) 

where -Fi+1/2 is the flux through the z i+1 / 2 interface. 

The results of closer coupling of the gravity to the hydrodynamics is shown in Figures 
2 and 6. In both figures, the results are from the simulation of an isothermal atmosphere, 
discussed in §3.1, with a pressure scale height at the base of 1036 cm and sound speed at 
the base of 4.4 x 10 8 cms -1 . In Figure 2 we show the effects of using LeVeque's method in a 
Godunov solver, and in Figure 6 we show the results of using the above- described modified 
states in PPM. In both cases, we can see a corresponding increase in accuracy of the solution. 

Note that neither the two figures, nor the two methods of dealing with the source-terms, 
can be directly compared, as the first figure demonstrates results with a very low order 
Godunov method, and the second from using higher-order PPM. The figures demonstrate 
only that our source-term correction to our hydrodynamic solver gives a similar improvement 
in results as other approaches in the literature work for other solvers. 



3. Initial Models 

To better understand the effects of the initialization, boundary conditions, and hydro- 
dynamics, we consider three different initial models. We choose the conditions of all of the 
initial models to roughly agree — a base density of a few x 10 6 g cm -3 and a base tempera- 
ture of ~ 10 7 K, with a solar-like composition. We use a constant gravitational acceleration, 
g = —1.8676 x 10 14 cm s -2 . These conditions are relevant to an accreted atmosphere on a 
1.4 solar mass, R = 10 6 cm neutron star. The accretion phase is computed using a one- 
dimensional implicit hydrodynamics code for many dynamical times, and the subsequent 
burning front propagation is then followed using a multidimensional hydrodynamics code. 

Once the initial model is set up, it will be mapped onto the computational grid. The 
PPM algorithm we use carries the cell averaged value of each variable in each zone. For 
an initial model that is not analytic, but rather a series of points and corresponding values, 
we will need to map onto the grid with the understanding that the values on the new grid 
are treated as zone averages. We will consider several different initial models of increasing 
complexity below. 
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3.1. Isothermal Ideal Gas EOS Atmosphere 

The simplest initial model one can imagine is an isothermal atmosphere in hydrostatic 
equilibrium, whose pressure is obtained via the ideal gas law, 

dP I I / \ p k B pT 2 

dF = -'^ ; P= A^ =PC - (33) 

where P is the pressure, g is the constant gravitational acceleration, p is the density, ks 
is the Boltzmann constant, T is the temperature, A is the average atomic mass, m u is the 
nucleon mass, and c s is the isothermal sound speed. Although this EOS is not relevant to 
a neutron star atmosphere, we will use it to examine some of the properties of hydrostatic 
equilibrium in a Godunov type code. 

Eq. (33) can be integrated analytically to yield an exponential atmosphere, 

p(z) = p -exp{-\g\z/c 2 s } = p -exp{-z/H} , (34) 

where po is the base density, and 

P(z) = p c 2 s exp{-z/H} = P exp{-z/H} . (35) 

The quantity H = c 2 /\g\, the scale height of the atmosphere, will play an important role in 
determining the resolution necessary to hold any of the initial models steady. If H < dz, 
where dz is the zone size, then the model will be very under-resolved, and we can expect to 
have great difficulty in maintaining HSE. 

To see this, consider modeling the exponential atmosphere with a parabola, similar 
to the way a second-order Godunov-type code will model the physical quantities. If the 
pressure at z = is Pq, and the left and right points are 5 on either side with cell-averaged 
pressures and (P) +1 (from Eq. 10), we can construct the quadratic that has the same 

cell-averaged quantities as the exponential; in particular, the derivative of pressure over the 
cell containing z = is 

'dP\ (P) +1 - (PU 

Tz/ Q = 25 • (36) 

The cell-averaged pressures on the left and right will be 

(P)-i = Poj (exp{-5/(2# )} - exp{-3cV(2tf )}) (37) 

and 

(P) +1 = Poj (exp{35/(2H)} - exp{5/(2H)}) , (38) 



-Ir- 



respectively. (This comes from integrating the exponential pressure, P(z) over the left and 
right cells, which range from (—35/2,-5/2) and (5/2,35/2) respectively.) With the above 
equation and Pq = pogH, and defining x = 5/H, the acceleration (in units of g) at z = is 



^(f) ="M cosh ©" cosh (!9) 



(39) 



which, to second order in x, is —1 — (5/24)x 2 . The first term of this acceleration is simply 
the gravitational acceleration, leaving an error in acceleration of magnitude ^x 2 . Per CFL 
timestep, this corresponds to a spurious velocity, expressed as a Mach number, of 

M ~ . (40) 

Thus, if one wishes to consider very subsonic flows M. ~ 10~ 3 moving through a stratified 
atmosphere, one needs at least 6 points per pressure scale height, or else one will get velocities 
of greater magnitude every timestep. Since the velocity errors are cumulative, one would 
realistically need much more than that. 



3.2. Realistic EOS Atmosphere 

The next complication we can imagine adding to our initial model is using a more 
realistic equation of state. This more generalized EOS prevents us from integrating the 
model analytically as above, but we can still compute the model on an arbitrary grid by 
simply differencing the equation of HSE (Eq. 13). To recover cell averaged quantities, we 
can subsample each zone and compute the average. 

We use a Helmholtz free energy, table-based EOS for a degenerate electron gas with 
perfect gas ions and radiation pressure included with FLASH (Timmes & Swesty 2000) for 
this model. The material at the base of this model is partially degenerate, and this EOS 
accurately describes the thermodynamics in our neutron star atmosphere. To complete this 
model, we again assume that our atmosphere is isothermal. The atmospheric structure is 
completely determined by the choice of temperature (we use 3 x 10 K), the composition (we 
assume 0.7 1 H and 0.3 4 He by mass), and the density at the base (we take 5 x 10 6 g cm -3 ). 
Because of the complexity of the EOS, despite the fact that this is an isothermal atmosphere, 
the scale height is not constant. 

The pressure and density profiles are computed as follows. At each point, the X± and 
T are taken as given, and the equation of hydrostatic balance, 

(41) 



13 



is numerically integrated outwards, using reconstructions to the data, as described below. 
This equation, plus the equation of state, is enough to determine (by iteration) both the 
pressure and density separately. 

Given the cell-averaged pressure and density data from the previous points and the 
functions in §2.1, we can construct functions for the pressure and density field up to and 
including the current cell. Given two such functions, one for density and one for pressure, 
we can use Eq. (41) to relate (p) +1 and (P) +v 

If we use a linear function to model the density, as in Eq. (9), we must use a quadratic 
one to model the pressure (because of the derivative), as in Eq. (10). Using these expressions 
in Eq. (41), and solving for (P) +1 and (p) +l , we get 

(P) + i ~ ( p )o = y «P> +1 + <P)o) , (42) 

which, with the EOS, is enough to determine (P) +1 and (p) +v We will refer to this expression 
as our first-order differencing. 

For more accuracy, we can use a quadratic to fit the density, and a cubic to fit the 
pressure. The cubic is given by Eq. (11). We then find 

(P) +1 ~ (P)o = ^ (5 <P> +1 + § (p) - (p>_J • (43) 

We will use both of these differencing schemes in §5 to assess how large of a difference the 
choice makes. 

A final difference scheme, using a quartic (Eq. 12) and a cubic, can be constructed: 

<^>+i - <^>o = ^ (9 <P> +1 + 19 (p) - 5 (pU + <pL 2 ) • (44) 

but this extra accuracy is not expected to be significant for a second-order accurate code, 
and numerical experiments confirm this. 

To generate our model, we take the base density and pressure and use either Eq. (42) 
or (43) to find the pressure and density in the next zone. We continue this procedure, 
moving outward from the base, until the density falls below a small density cutoff we impose 
(1CT 5 g cm -3 ). Figure 7 shows a plot of this model. We notice that as the density changes, 
the scale height also changes. At the top of the atmosphere, the small scale height is reflected 
in the steepness of the density just before our low-density cutoff. The sound speed at the 
base of the atmosphere is ~ 4 x 10 8 cm s _1 , giving a dynamical timescale for the atmosphere 
of ~ 5 ps. Figure 8 shows the relative error in the density of the model when using the 
first-order vs. second-order differencing. 
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3.3. Model from a 1-d Stellar Evolution Code 

The final model we consider is one that comes directly from the Kepler 1-d stellar evo- 
lution code (Weaver, Zimmerman, & Woosley 1978). This model was generated by accreting 
H/He onto the surface of a neutron star, gradually building up a fuel layer 17 meters thick. 
Nuclear burning, a complex equation of state, and convection were included in this calcula- 
tion. The evolution was stopped shortly before runaway, and the velocities in the atmosphere 
are small, but non-zero. Additionally, this code solved the hydrodynamics equations in a 
Lagrangian formulation, requiring a translation from mass-based zones to an Eulerian grid. 
The number of grid points in the Kepler model is much smaller than the number of points 
on our Eulerian grid. 

We note that in contrast to the other models described above, the composition in this 
model is not uniform, but abruptly changes at the transition from the accreted fuel layer 
to the underlying neutron star, which is comprised of iron. This in turn creates a density 
jump, since the ionic component of the pressure scales as l/A. 

To import this model onto the FLASH grid, we use the following prescription. First, the 
model is re-gridded from the original Lagrangian mesh onto a uniform, one-dimensional Eu- 
lerian grid, whose resolution is equal to the finest spatial resolution on the multidimensional 
hydrodynamic grid. Once interpolated, we renormalized the abundances. Varying precision 
in the initial model data file compounded by interpolation errors may result in a set of abun- 
dances in a zone that do not sum to unity; we divide by the actual numerical sum of the 
abundances to enforce the constraint of mass fractions summing to one. Next we take the 
temperature, density, and composition as given, and update the remaining thermodynamic 
variables in the zones with the FLASH EOS. In practice, the change in the thermodynamic 
variables caused by the differences in the EOS is small, < 1%. 

We now want to restore hydrostatic equilibrium to this model, after adjusting the ther- 
modynamic variables. The velocities are set to zero in every zone, since we are not presently 
concerned with how to map a one-dimensional velocity field into multiple dimensions. If the 
model was really from a slowly simmering phase, this is not a bad approximation; in the 
model used in these calculations, for instance, the maximum velocity was ~ 10 -1 cm s -1 . To 
restore this model into hydrostatic equilibrium, we must pick a point in the model whose p, 
T, and Aj will remain fixed, and integrate outward from there. We show two choices here, 
(i), taking the bottom-most zone as the starting point, and (ii), taking the point just above 
the composition interface as the starting point. 

The differencing is performed in the same manner as described in the previous section. 
We experiment with both differencing schemes, Eqs. (42) and (43). This differencing is 
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continued throughout the entire model. When integrating toward the top of the atmosphere, 
we stop putting the model into hydrostatic equilibrium once the density becomes so small 
as to be no longer dynamically important. We use a cutoff value for the density of 1 x 
10~ 5 g cm" 3 , and affectionately refer to the material above this as the "fluff". This cutoff is 
needed since we cannot continue the HSE profile to arbitrary heights, as the densities would 
quickly underflow. 

We note that instead of adjusting the density along with the pressure, one could adjust 
the temperature. For cases where the atmosphere is degenerate (like our present case), this 
can be problematic due to the insensitivity of the pressure on the temperature. We do not 
explore this approach in the present paper. 

Figure 9 shows the results of applying this procedure to our Kepler initial model. The 
two panels differ in the choice of the reference point used when differencing the model 
into HSE. We note that in the top panel, where we chose the base of the model to begin 
the integration, the errors compound greatly as we integrate outward, especially where we 
have to integrate through a material discontinuity, which is poorly modeled by a low-order 
polynomial. At the base of the fuel layer, there is a significant difference in the density. The 
agreement is much more uniform in the lower panel, since the reference point (the base of the 
fuel layer) is near the center of the model. We use the base of the fuel layer as the reference 
point in all Kepler-model simulations presented here. 

4. Boundary Conditions 

The stability of the model atmosphere can be very sensitive to the choice of boundary 
conditions. In these highly-subsonic simulations, the hydrodynamic equations are essentially 
elliptic, so that boundary conditions matter as much as the initial conditions inside the 
computational domain. We investigate many different boundary conditions at the lower 
boundary, and two different boundary conditions at the upper boundary in this study. 

In FLASH, like most finite-volume hydrodynamics codes, the boundary conditions are 
implemented in fictitious zones outside of the physical domain called guardcells or ghostcells. 
In order to allow for refinement and parallelization of the code, the computational domain 
is broken into multiple sub-domains, or blocks. Each block is surrounded by a perimeter of 
guardcells that hold the data from neighboring blocks, or, if at a physical boundary of the 
computational domain, are filled with the proper boundary condition. 

The problem of generating good hydrostatic boundary conditions is closely related to 
that of finding good 'outflow' or non-reflecting boundary conditions (Thompson 1987), which 
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remains an area of research. In both situations, to create a desired flow condition, one has to 
set up corresponding physical fluid conditions, essentially inverting the Riemann problem. 
Further compounding the difficulty is that the problem must be solved in a way consistent 
with one's hydrodynamic solver. Below are several boundary conditions, both commonly- 
used in the astrophysical literature and novel, which are approximations of solutions to the 
inverse problem. 

4.1. Lower Boundary 

The lower boundary must support the weight of the atmosphere above it, while still al- 
lowing for dynamics. We consider two classes of boundary conditions — a standard reflecting 
boundary, and a hydrostatic boundary, which provides pressure support to the material above 
while still allowing flow through the boundary. 

4-1.1. Reflecting 

The simplest lower boundary we can use that will support the weight of the material 
above it is a reflecting boundary. This is one of the most commonly used boundary conditions 
when evolving a hydrostatic atmosphere. Absent any gravity, this boundary condition simply 
reverses the sign of the normal velocities in the boundary region, and gives all other variables 
a zero- gradient. The result is that there is no flow through the boundary. Any wave that 
hits the boundary is reflected back into the interior. This boundary condition is effective at 
restricting flow through the boundary, but it will not allow sound waves to leave the grid as 
the initial model relaxes. 

In the presence of gravity, the traditional reflecting velocity boundary condition of 'flip- 
ping' the velocity in the direction transverse to the boundary, 

v z (z - z ) = -v z (z - z) , (45) 

(where zq is the location of the physical boundary) will not work, because there is an accel- 
eration term due to gravity. Better is 

v z {z - zq) = -v z (z -z) + g5t , (46) 

which takes into account the acceleration performed at the hydro step when integrating 
the equations. In a finite volume-code it is easy, in addition, to enforce the desired 'no- 
penetration' condition exactly by setting the flux across the interface defining the boundary 
to exactly zero. This is done in the reflecting boundary condition results below. 
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4-1.2. Hydrostatic 

The reflecting boundary is artificial in the sense that it does not allow waves to flow 
off the grid. An alternative is to use a boundary condition that understands hydrostatic 
equilibrium. The basic strategy is to provide pressure support to the material above the 
boundary by solving the equation of hydrostatic equilibrium in the boundary region. This 
is done by simple differencing as described above. When filling the guardcells, an additional 
constraint on the material is needed, since the HSE equation does not uniquely determine the 
boundary states given the states above. We assume that the composition and either density, 
temperature, or entropy are constant. This is not a complete set of possible constraints, 
but is enough to illustrate the effect this can have on the evolution of an atmosphere. Any 
constraint should be motivated by the physics of the model under study. Figure 10 shows 
the density, temperature, and pressure at the base of the Kepler model (including guardcells) 
for the different choices of constant variable. We see that they all yield roughly the same 
pressure profile, since the pressure at any point in the atmosphere is determined by the 
weight of the material above it: 

P = J pgdz = ga , (47) 

where a is the column density. Because of the degeneracy of the gas, we see that the 
temperature must rise dramatically in the guardcells in the constant density case in order 
to provide the needed pressure support. 

It is important to note that the best boundary condition may be problem dependent. 
Differences in the EOS or the physics of the atmosphere will affect the choice of the constraint. 



constant p, first order 

The easiest hydrostatic boundary to implement uses the simple first-order differencing of the 
hydrostatic equilibrium equation (Eq. 42). In §3.2, we derived this to find the pressure in 
the zone immediately above the current zone. At the lower boundary, we need to find the 
pressure in the zone below the last interior zone, given the pressure and density in that zone. 
Thus, our difference equation becomes: 

( P )o - = - y {{p)o + {p)+i) ■ (48) 

The right hand side requires the density in the th zone. In this case, we assume that the 
density is constant (zero-gradient) in the guardcells, initializing all of them with the value of 
the density in the first interior zone. Therefore in order to satisfy the EOS, the temperature 
will rise in the guardcells. 
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All other variables are also given a zero-gradient, except the velocity. We use three 
different methods for dealing with the velocity. The first is to give it a zero gradient, so the 
velocity in the guardcells is constant, and equal in magnitude and direction to the velocity in 
the last interior zone. We call this the "outflow" condition. The next method is to perform 
the outflow method only if the velocity is leaving the grid. If the sign of the velocity in 
the first interior zone is positive, indicating that material is flowing onto the grid, we set 
the velocity in the guardcells to 0. We call this the "diode" boundary condition. The final 
method is to simply reflect the velocity, as indicated by Eq. (45). 



constant T, first order 

A slight variation on the above is to give the temperature a zero-gradient in the guardcells, 
initializing it to the value in the first interior zone. Then the density will need to increase (as 
given by the EOS) in order to give the pressure demanded by Eq (48). We also experiment 
with the three different velocity methods described above. 



constant entropy, first order 

A more realistic boundary condition for a stellar model is to make the temperature/density 
profile isentropic and in HSE. For an atmosphere in hydrostatic equilibrium the adiabatic 
temperature gradient is 

(t) =— ' (49) 

\ dz Jad °P 

(see, for example, Lantz & Fan 1999) where 

d hip 



<91nT 



(50) 



and c p is the specific heat at constant pressure. Together with the equation of hydrostatic 
equilibrium, this specifies the conditions in the guardcells. We still have to assume a function 
for the velocities, and we use the same choices as above. 

To implement this boundary condition, we use a first order differencing for the HSE 
equation (Eq. 48), and a simple first order differencing of Eq. (49). These two expressions 
are solved simultaneously along with the EOS and iterated until we obtain convergence of 
the pressure and density. 
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constant p, second order 

We can use a variation of Eq. (43) to perform the differencing in the guardcells. This is 
a higher order differencing than that used in the above boundary conditions, and should 
give more accurate results. Again, we are interested in using the information in the zones 
just above the lower edge of the computational domain to fill the guardcells below it, so the 
differencing is in the opposite direction as before: 

(^)o " <^>+i = ~ (5 (P) + 8 (P) +1 - (p) +2 ) • (51) 

Again we give all variables except velocity a zero-gradient, and retain our three choices for 
dealing with the velocity as described above. In this case we keep the density constant in 
the guardcells and the pressure given by Eq (51) as a constraint on our EOS to give the 
temperature in the guardcells. Thus the temperature will rise with depth in our boundary. 

constant T, second order 

This case is the same as above, but we take the temperature as constant in the boundary 
zones and solve Eq (51) in tandem with the EOS to find both (P) and (p) simultaneously. 
Again, we allow for all three different velocity choices. 

constant entropy, second order 

This case is another implementation of the constant entropy boundary condition defined 
above, but we use the second order differencing for the HSE equation (Eq. 51), and keep 
the first order differencing for the temperature expression. The pressure and density are the 
two variables that are involved in the dynamics, so this is why these are treated as second 
order. 



4.2. Upper Boundary 

The choice of the upper boundary is also important. Further out in the hydrostatic 
envelope, eventually the density will reach a very small value, such that it can no longer be 
represented in IEEE floating point arithmetic without underflowing. Figure 7 illustrates this 
effect nicely. Above 2000 cm, the density drops off very rapidly. If we were to continue to 
follow the density down to arbitrarily small values, the change in the height of the atmosphere 
would be insignificant, to the point where it would be less than a computational zone. There 
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are two major paths one may take to get around this. First, one may only put a portion of 
the atmosphere on the grid, leaving off the top few scale heights. In this case, a hydrostatic 
boundary condition is required at the top, preferably one that makes the same assumptions 
as the lower boundary. The difficulty with this method is that, if in the dynamics that follow 
in a simulation the envelope heats up, the scale height will get larger, and more and more 
of the envelope will be lost through the top boundary. This can be overcome partially by 
employing an expanding grid, but that case is not considered here. Once material leaves the 
top of the grid, it is lost forever, and would be unable to further contribute to the dynamics 
of the simulation. 

An alternative is to follow the hydrostatic structure down to a cutoff value of density 
(something large enough still so it will not underflow), and then apply a uniform density 
above this point in the atmosphere. This creates the 'fluff' described in §3.3. The top 
boundary can then be a zero gradient, or impose some inflow characteristic of accretion. Since 
we do not place this material in hydrostatic equilibrium, it will fall under the influence of 
gravity. However, its mass is so small that it will have very little dynamical impact. Leaving 
a buffer of the low-density material above the hydrostatic envelope allows the envelope to 
expand in the computational domain. In particular, expansion that is non-uniform can 
more easily be handled in this case than with an expanding Eulerian grid. This choice of 
boundary condition allows the exploration of surface features, such as waves, that are created 
as localized heating causes some regions of the envelope to expand before others. We choose 
this avenue for most of the simulations presented in this paper, and study the effect of this 
choice in 55.1.2. 



5. Results 

We ran simulations using both the isothermal, realistic EOS atmosphere, and the Ke- 
pler generated initial atmosphere. Both initial models were run many times, varying the 
resolution, grid type, boundary conditions, and hydrodynamics. Unless otherwise noted, all 
simulations use the standard PPM algorithm. Table 1 summarizes the different parameters 
we explore. All calculations were performed in one dimension, for computational efficiency. 
We examine the effect of each of these parameters in turn below. 
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5.1. Effect of Boundary Conditions 

5.1.1. Lower BC 

As discussed in §4.1, there are a large number of variations of the hydrostatic boundary 
condition. We consider the effect of the boundary conditions first, since, as we will see, 
the choice of boundary condition has a large effect. To examine the influence of all these 
differences, we ran each of our initial models on a uniform grid with the standard PPM 
algorithm, varying the parameters in the boundary conditions. All of these calculations 
have a zone width of 2.4 cm and were run for 250 fis (~ 20 dynamical times). The results 
are summarized in Figures 11 and 12. 

A total of 10 cases are examined here for both the Kepler initial model and the isothermal 
atmosphere initial model, showing both the density and velocity profiles at the end of the 
calculation (250 fis), the total momentum as a function of time, and the kinetic energy 
as a function of time. These latter two plots allow us to understand how the velocities 
are evolving with time — whether they are relaxing down to a quiet state (as we hope), or 
increasing without bound. Reflecting and hydrostatic boundaries were used at the lower 
boundary, and a fluff condition was used at the top. We note that this momentum includes 
that of the fluff, but since its density is so low, we expect it to only make a minor contribution. 
Despite our best efforts at initializing the grid with our model in a manner consistent with 
our choice of hydrodynamics algorithm, some velocities are quickly generated as we evolve 
our models in time. All simulations were initialized using the second order differencing 
and the second order variants of the different boundary conditions were used. Ideally, the 
velocities would be small in magnitude and diminish quickly as the model relaxes. We will 
look at the effect of the initialization method in a latter section. 

We note a number of things right away — the hydrostatic boundary conditions that keep 
the temperature or entropy constant do a far better job than those that keep the density 
constant. In Figures 12d-e, we note that the Kepler atmosphere is falling off the bottom 
of the grid and large negative velocities are dominating in the atmosphere. In the case of 
hydrostatic with constant density, but reflecting velocities (Figure 12f), we note that the 
atmosphere is held in the box, due to the reflection of the velocities. We refer back to Figure 
10 — the extreme rise in temperature in the guardcells needed to support the weight of the 
material above the boundaries is the likely culprit here. This situation is much more severe 
for the isothermal atmosphere case (Figure lld-e), as the model quickly (< 50 /is) falls 
through the bottom of the domain. 

The hydrostatic conditions that use a simple zero-gradient /outflow condition on the 
velocities for the constant temperature case (Figures 11a and 12a) show the momentum 
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monotonically increasing with time, while the constant entropy case (Figures llg and 12g) 
show the momentum monotonically decreasing with time. The velocity is very flat in the 
hydrostatic atmosphere, with the exception of a 'hiccup' at the composition interface of the 
Kepler model. However, this velocities steadily increases in magnitude with time, making 
this condition ill-suited for long simulations. 

Any condition (hydrostatic or not) that reflects the velocities at the lower boundary 
(Figures llc,f,i-j and 12c,f,i-j) shows a ringing which may be observed in the momentum and 
energy plots, with a period about equal to the dynamical timescale for the envelope. This 
is due to the reflection of the velocities, which, with the hydrostatic state, diminishes the 
waves that can penetrate through the boundary. 

The diode constraint coupled with the constant temperature hydrostatic boundary con- 
dition also shows a ringing, with an amplitude (see Figures lib and 12b) that is on the order 
of, or less than that of the reflecting boundaries (compare to Figures 11c and 12c). It also 
appears to be decreasing in amplitude with time. After an initial transient, the diode bound- 
ary conditions give the lowest velocities in the envelope of any of the boundary conditions 
tested. 

We can also look at the effect of the order of the HSE differencing in the boundary 
condition on the velocities in the envelope. In all cases, the velocities and magnitude of 
the momentum are smaller with the 2nd order differencing (holding the other parameters 
constant). 

The boundary condition that leads to the quietest velocity field, after the initial tran- 
sients die down, is the hydrostatic boundary with constant temperature and diode velocities. 
For the Kepler model, this boundary condition gives a momentum after 250 fis that is 2 or- 
ders of magnitude smaller than the reflecting boundary. We stress again that this finding 
may be problem dependent, and tests should be done on any new problems to see if this 
remains the case. 



5.1.2. Upper BC 

Figure 13 shows the density and velocity as a function of height at several times for the 
isothermal atmosphere initial model, with both the fluff upper boundary and a hydrostatic 
boundary at the top of the domain. The hydrostatic boundary was first order, constant 
temperature, with a diode condition on the velocity. The isothermal, complex EOS initial 
model was used for these runs. The two runs were done with the same spatial resolution, on 
a uniform mesh. In the fluff case, the domain extended to 4000 cm, but we only show the 
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lower 2000 cm, to match the domain used in the hydrostatic upper boundary case. Figure 
14 shows the same comparison for the Kepler initial model, again run at identical spatial 
resolutions, with domain sizes of 3000 cm (without fluff) and 6000 cm (with fluff). 

In both figures we see that the fluff boundary at the top of the domain is at least as 
effective as a hydrostatic boundary at the top in maintaining a hydrostatic atmosphere. In 
fact, the velocities generated throughout the atmosphere with the fluff boundary condition 
are generally smaller in magnitude than those with the upper hydrostatic boundary. Coupled 
with the benefit that the envelope is able to expand on our grid with out losing dynamically 
important mass through the top boundary, the fluff condition at the upper boundary is the 
optimal choice. 

5.2. Effect of Initialization Method 

All of the boundary condition comparisons shown above were done with second order 
differencing in the initialization. We can also use the first order differencing (Eq. 42) to put 
the initial model into HSE on our grid. This does not have nearly as large an effect as the 
different boundary conditions do. Figure 15 shows the Kepler initial model with both first 
and second order differencing done at initialization. Both runs used the hydrostatic lower 
boundary with constant temperature and diode velocities, and the standard PPM algorithm. 

Looking at the plot of momentum verses time, we notice that at early times, the first 
order case generates higher positive velocities initially than the second order the 
model relaxes to the grid. The two cases quickly settle down to roughly the same state after 
about 50 /is, and by the end of the calculation (250 fis), the plot of velocity as a function 
of height does not show much difference between the two runs. This is to be expected; 
absent any driving terms from poor boundary conditions or energy inputs, the simulation 
will eventually 'find' hydrostatic equilibrium on its own. However, a good choice of initial 
model will strongly reduce the initial transients. This model is more complicated than the 
simple isothermal atmosphere because of the composition interface. This leads to the more 
severe velocities. We note that this comparison was done with the best choice boundary 
conditions we had available. If the boundary conditions are poor, or the problem is under 
resolved, then spurious velocities will continue to be driven in the simulation every timestep. 
Everything else, no matter how bad it is, will eventually even out as the Euler equations do 
their work and find HSE. The order of the initialization method should match the order of 
the boundary condition in to avoid transients caused by jumps at the boundaries. 

In some problems, the initial transient may not be a problem. In others, however, they 
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may interfere with a phenomenon one is trying to measure, or — if the transients cause 
spurious burning, or cause material unphysically to leave your computational domain - 
they may alter or change entirely the long-term evolution of the system. 

5.3. Effect of Resolution 

Figure 16 shows the kinetic energy verses time for an isothermal atmosphere with 2nd- 
order differencing and 2nd-order constant-temperature boundary conditions. This figure 
shows only the short term evolution. Our estimate in §3.1 showed the importance of re- 
solving the scale height in minimizing any spurious velocities. We therefore expect to see 
improvement in the ambient velocities as we increase the spatial resolution. As we might 
expect, our model shows second-order convergence (see Figure 17) in resolution. 

We can combine the effects of higher-order initialization and resolution by looking at the 
convergence behavior of the velocities in the isothermal atmosphere with different initializa- 
tion schemes. This is also shown in Figure 17. We see that with the operator-split gravity, 
we still only get 2nd order convergence even with the higher-order spatial discretization for 
the initial model, however the spurious kinetic energy is an order of magnitude lower. 

Figure 18 shows the long term evolution of the Kepler initial model for four different 
spatial resolutions. The '6 level' curve corresponds to that shown in Figure 12. All four sim- 
ulations were run with the first order initialization and the constant temperature hydrostatic 
boundary conditions. We see that as the resolution is increased (each level is a doubling in 
resolution), the momentum and kinetic energy decreases dramatically, as we would expect. 

5.4. Effect of Input States to Riemann Problem 

We performed, with our modified-states PPM, the same numerical experiment as in §5.3, 
which already had fairly low velocities. The results are shown in Figure 19. By modifying 
the input states, the numerically-induced kinetic energy drops by three orders of magnitude 
over those shown in Figure 16, and converges to third-order, rather than to second-order, as 
shown in Figure 20. The more direct coupling of the gravity into the hydrodynamical solver 
increases both the accuracy of the solution and the convergence properties. 

This improvement will only become apparent when using a good initialization method 
and boundary condition, otherwise, any improvement it may yield will be swamped by noise 
from these other sources. We used the second order constant temperature/diode velocity 
boundaries for these simulations to minimize any boundary effects. 
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6. Conclusions 

We have studied the effects of the initialization process, boundary conditions, resolution, 
and the hydrodynamic algorithm on the process of mapping an initial model onto a Eulerian 
grid in a Godunov-type code. We saw that depending on how one goes about the process, the 
result can be a nicely relaxed, quiet mapping, or one dominated by high velocities, swamping 
out any physical processes that one may be interested in studying. 

Boundary conditions have the greatest impact on maintaining the stability of the hy- 
drostatic atmosphere. The standard reflecting boundary is poorly suited for maintaining a 
hydrostatic envelope — the inability of minor pressure disturbances to escape off the bottom 
of the grid as the model relaxes results in high velocities. The actual choice of boundary 
condition should reflect the physics of the atmosphere being supported. If the initial model 
is isothermal, then an isothermal/HSE boundary condition provides with best match. The 
constant entropy boundary may be a better match for an atmosphere that is generated by a 
stellar evolution code. For any model, a hydrostatic boundary condition, using differencing 
that matches that of the initialization, with a secondary constraint that matches the physics 
of the model atmosphere is the best solution. Assuming that the density is constant in 
the boundary proved to be the worst assumption. The degeneracy of the EOS puts strong 
demands on the temperature profile to counteract the weight of the atmosphere. This as- 
sumption may fare better with a gamma-law EOS, but this was not considered in the present 
paper. 

We also demonstrated that, once an appropriate boundary condition is chosen, the 
treatment of source terms in the hydrodynamic solver can have an impact on the stability 
of the atmosphere. Our modified-states PPM was effective in reducing the magnitude of 
the spurious velocities generated as our model relaxed onto our grid. This method, or that 
proposed by LeVeque (1998b), can be adopted to most Godunov type codes to increase the 
accuracy of the hydrodynamics in the presence of gravity. 

As expected, the resolution is important in reducing the errors in maintaining a hydro- 
static atmosphere. The resolution should be chosen to be fine enough to resolve the scale 
height of the atmosphere well, and to keep ambient velocities smaller in magnitude than 
whatever physical processes under study would yield. In an AMR calculation, the coarsest 
resolution used in the simulation must be this fine, even in regions of the atmosphere with 
no features. 

Boundary conditions, source terms, and poor resolution can all be sources of spurious 
velocities throughout a calculation. By contrast, initialization methods can cause at most 
a transient while the simulation settles into hydrostatic equilibrium. This transient may 
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take a while to settle, however, or may have dynamical consequences later in the simulation, 
so care must be taken. The initialization methods described in this paper produce profiles 
which generate very small transients. 

We note that the methods that we describe in this paper can be applied to any initial 
hydrostatic atmosphere. As a final example, we show the havoc a poorly initialized model can 
wreak in a simulation. Figure 21 represents two-dimensional simulations evolved from a one 
dimensional model provided by S. A. Glasner (2002, private communication) of a classical 
nova precursor; it is a model from the same simulation that produced the Id model used in 
Glasner, Livne, & Truran (1997); Kercek, Hillebrandt, Truran (1998), but at an earlier time, 
before convection begins. Thus, one could hope to examine the multidimensional onset of 
convection in these nova precursors (Dursi et al. 2002). The convection is driven by nuclear 
reactions 'simmering' near the interface between the C/O white dwarf and the accreted layer 
of stellar material. Since the vertical scale of the simulation is significant compared to the 
radius of the white dwarf, plane-parallel inverse-square gravity is used rather than constant 
gravity. 

The model is perturbed in the highest-temperature region of the accreted atmosphere 
with a 10% temperature increase at time t — 0. The white contours in the figures enclose 
the region in the simulation with a temperature greater than this perturbed temperature. 
The black contour marks the white dwarf /accreted material interface. 

Without taking any of the precautions outlined in this paper, an unphysical 'settling' 
occurs, caused by poor boundary conditions and differences in the EOS. This generates 
large velocities («~5x 10 6 cms -1 ) and compressional heating, as shown in the figure. The 
heating, combined with unphysical mixing across the interface caused by the large velocities, 
then cause a completely spurious layer of increased burning, which then dominates the long- 
term evolution of the simulation. By contrast, the simulation using constant-temperature 
second-order boundary conditions, second-order initialization, and the modified-states PPM 
shows the beginning of the formation of convective rolls. Both simulations were done on a 
1920 x 640 uniform mesh, with a computational domain of 1 x 10 7 cm x 3 x 10 7 cm; Figure 
21 shows only the domain near the interface of the white dwarf and the accreted material. 
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Appendix: Non-Constant Gravity 

In this paper, we have considered only the case of constant gravity. It is relatively 
straightforward to extend this to the case of fixed but non-spatially-constant gravitational 
acceleration. 

Consider for instance Eq. (42). This was derived by using a linear reconstruction for 
the {p) i data and putting into the equation of hydrostatic balance (Eq. 41). If instead, the 
quantity (pg) is reconstructed, and the same procedure applied, one ends with 

(P} +1 -(P) = 6 -{(pg) +1 + (pg} ) . (52) 

However, to iterate on the equation of state and to calculate a local average density, one now 
must deconvolve (pg) 1 to find (p) 1 . There are a number of ways one could do this. One way 
is to use a cell-averaging algorithm consistent with PPM: 

{f) = \[fL + ±fc + f R ] , (53) 

where the subscripts (L, C, R) refer to the values of the reconstructed function at the left cell 
interface, cell center, and right cell interface, respectively. Presumably, whatever our fixed 
gravitational acceleration, we can evaluate it numerically at any given point, so that we can 
then use 

(54) 



(P9)i = I 



with the computed values of g and evaluate the pointwise values of p using the (here, linear) 
reconstruction of density to solve for the local density average: 

/ \ _ (P)o (9r ~ 9l) + 12 (pg) 1 , . 

{P)1 ~ 8 9c + 9l + 39r " ( } 

One can repeat this procedure from Eq. (54) using the quadratic reconstruction to get 

, . (p) ggg - 2g C ~ 5g L ) + {-2g R + g c + g L ) + 36(pg) 1 

{Ph ~ 23g c + 2g L + llg R ' 1 j 

Thus, one can use the first- or second-order differencing in both the initialization and 
the boundary condition to find (pg) ± , and then the above equations to find the implied value 
of (p) v and iterate on the equation of state as before. 
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Table 1: Summary of the different parameters used in the grid of calculations 
property value 
initial model isothermal 

kepler model 
initialization method first order differencing 

second order differencing 
hydrodynamics ppm 

ppm/modified-states 
levels of refinement 5 

6 

boundary conditions reflect 

constant p, first order 
constant T, first order 
constant s, first order 
constant p, second order 
constant T, second order 
constant s, second order 
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Fig. 1. — A function po(z) is computed, which reconstructs the data for p within zone 0. 
Shown here is a cubic, as in Eq. (11); it averages over Zi = (—25, —5,0,5) to the previous 
data values ((p)_ 2 > (p)-i > (p)o > (p)i)- Reconstructions using lower orders of polynomials 
need fewer data points. 
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Fig. 2. — Velocity profile, as a function of vertical position, of the isothermal atmosphere at 
time t = 2/is (approximately 1/2 of a dynamical time) when evolved with 3072 points with 
Godunov's method (solid) and Godunov's method corrected by LeVeque's method (LeVeque 
1998b) (dashed). The maximum velocity in the corrected case, which cannot be seen on the 
scale of the graph, is 975 cms -1 . The atmosphere's pressure scale height ranges from 1036 cm 
at the base to 48 cm at the top boundary; the sound speed ranges from 4.4 x 10 8 cms -1 to 
9.5 x 10 7 cms -1 . 
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Fig. 3. — The PPM reconstruction process for a simple exponential atmosphere, P(z) = 
Poexp{—z/H}, where H is the scale height. We choose c s = Pq = 1, and g = —4. The 
left column shows the PPM reconstruction on the pressure and the right column shows the 
reconstruction on pg, which is used to derive the wave-generating pressure. The top row is 
the analytic plots of the pressure and pg. The second row is the zone average values of each 
quantity (horizontal dashed lines) and the interface values (diamonds), as determined using 

the PPM internnlants The 7,nne enVes are marked with vertical Hotted lines The bottom 
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Fig. 4. — The wave generating pressure (solid lines) for the atmosphere shown in Figure 3 
as computed using Eqs. (27) and (28). The zone averages of the original pressure are also 
shown (dashed lines). 
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Fig. 5. — Another, 'close-up' view of the process of modifying the pressure used in calculating 
the states to the Riemann problem. Shown in (a) is an exponential atmosphere, P = 
Poexp{—z/H}, with Po = c s = 1, g — —4, and 5z = 1/15, with both the PPM reconstruction 
as in Eq. 16 (thick line) and cell-averaged values (dashed). In (b) is shown only the averages 
and reconstructions — as given by Eqs. 25 and 26 — of the cell to the left and right of 
the interface at z — 5z. Plotted in (c) is the wave-generating pressure on either side of the 
interface, as in Eqs. 27 and 28. As can be seen, this modified pressure field for calculating the 
left and right states is essentially flat to the (quadratic) order of the original reconstruction; 
here, the third-order term is dominant. In (d) is the wave-generating pressure and the full 
pressure plotted together. 
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Fig. 6. — Same as Figure 2, but using PPM (solid) and our Modified-States PPM (dashed). 
The maximum velocity in the corrected case is 0.35 cms" 1 . The difference in magnitude 
between the velocities in the two figures is due to the difference in the accuracy of the 
hydrodynamic solvers (Godunov vs. PPM), not the result of the method used to deal with 
the source-term. 
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Fig. 7. — The density profile of an isothermal, uniform composition initial model with a 
complex EOS. Above 2000 cm, the density rapidly falls. 
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Fig. 8. — Relative difference in the density structure for the first-order vs. second order 
differencing. The error is greatest just before the low density cutoff. 
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Fig. 9. — Results of taking a 1-d initial model from the Kepler implicit stellar hydrodynamics 
code and putting it into HSE with the new EOS. We note the jump in density in this model, 
owing to the abrupt change in composition from the underlying neutron star to the accreted 
fuel layer. In the top panel, we took the base of the initial model as the reference point. In 
the lower panel, we took the base of the fuel layer in the initial model as the reference point. 
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Fig. 10. — Comparison of the three different hydrostatic boundary conditions, showing con- 
stant density (left), constant temperature (center), and constant entropy (right). Pressure, 
density, and temperature are plotted, scaled to their base values (x=0). Negative values of 
x are the guardcells. The temperature in the constant density case grows very rapidly in the 
guardcells, due to the degeneracy of the gas. 
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y 2nd order hydrostatic with constant temp and outflow velocities 
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b) 2nd order hydrostatic with constant temp and diode velocities 
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c) 2nd order hydrosiatic with constant temp and reflect velocities 




Fig. 11. — Effect of the boundary conditions on the velocity for the isothermal initial model. 
Each simulation was run for 250 fxs with a spatial resolution of 2.4 cm. The left panel shows 
the density and velocity as a function of height at 250 fis. The middle panel shows the total 
momentum as a function of time, and the right panel shows the kinetic energy as a function 
of time. 
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d) 2nd order hydrostatic with constan: dens and outflow velocities 
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e) 2nd order hydrostatic with constant dens and diode velocities 




f) 2nd order hydrostatic with constant dens and reflect velocities 
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11. — cont. 
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g) 2nd order hydrostatic with constant entr and outflow velocities 




h) 2nd order hydrostatic with constant en:r and diode velocities 





11. — cont. 




Fig. 11. — cont. 
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Fig. 12. — Effect of the boundary conditions on the velocity for the Kepler initial model. For 
each choice of boundary condition, three plots are shown: a snapshot in time of the velocity 
structure vs. height (taken at 250 fis); a plot of the total momentum vs. time, and a plot of 
the total kinetic energy vs. time. 
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d) 2nd order hydrostatic with constan: dens and outflow velocities 




e) 2nd order hydrostatic with constant dens and diode velocities 




f) 2nd order hydrostatic with constant dens and reflect velocities 
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g) 2nd order hydrostatic with constant entr and outflow velocities 
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h) 2nd order hydrostatic with constant en:r and diode velocities 
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Fig. 12. — cont. 
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120.003 fis " 240.001 fis 




Fig. 13. — A comparison of the velocities resulting when the top boundary is hydrostatic 
(dashed line) verses a buffer of low density 'fluff' (dotted line) for the isothermal atmosphere 
initial model. The four panels show the profiles at 120, 240, 360, and 480 /is. We see that 
the 'fluff' case yields lower velocities throughout the atmosphere. 
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Fig. 14. — A comparison of the velocities resulting when the top boundary is hydrostatic 
(dashed line) verses a buffer of low density 'fluff' (dotted line) for the Kepler initial model. 
As in the simple isothermal model, the 'fluff' boundary condition yields lower velocities 
throughout the atmosphere. 
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b) 2nd order initialization 




Fig. 15. — Comparison of the order of the HSE differencing used at the initialization stage 
for the Kepler initial model. Both runs used a hydrostatic lower boundary with constant 
temperature and diode velocity condition. Note that the higher-order initialization has a 
much (~ x8) quieter initial transient, but that medium- and long-term evolution is identical. 
This is as expected unless there is physics (eg., burning) in the simulation whereby the initial 
transients can feed back into the long-term evolution. 
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Fig. 16. — Plots of kinetic energy for the isothermal atmosphere described in §3.2. Since this 
model should be hydrostatic, the 'correct' KE is zero. The second-order differencing scheme 
was used for initialization, and second-order constant-temperature boundary conditions were 
used with diode velocity conditions. For scale, the mass in the simulation domain is ~ 
3 x 10 9 g, so that a KE of 10 15 erg corresponds to a density-weighted RMS velocity of 
6 x 10 2 cm/s, while the sound speed of the medium is from 1 — 4 x 10 8 cm/s. 
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Fig. 17. — The maximum kinetic energy in Figure 16 (dashed), and with the corresponding 
simulations run with first order BCs and initial differencing (solid), both with PPM. We see 
in both cases the kinetic energy converging as the 4 th power of resolution, meaning that the 
velocities are converging as the 2 nd power. 
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a) 5 levels of refinement 
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Fig. 18. — Effect of the spatial resolution on the velocity for the Kepler initial model. Four 
different resolutions are shown, a) 4.8 cm, b) 2.4 cm, c) 1.2 cm, and d) 0.6 cm. 
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c) 7 levels of refinement 




- S -2x10 



-6x10'* 

1000 2000 3000 4000 5000 5.0x1 0" 5 I.OxlO -4 1 .5* 1 0~ 4 2. Ox 1 0" 4 2.5* 

x [cm] tirre [V 




5.0x10 J 1.0x10^ 1.5x1 



2.0x10 2.5x10 



d) 8 levels of refinement 
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Fig. 18. — cont. 
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Fig. 19. — Same as Figure 16, but with the simulations run with the modified PPM states 
described in 82.2. 
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Fig. 20. — Same as Figure 17, but with the simulations run with the modified PPM states 
described in §2.2 added. We see the second-order convergence turning into third-order in 
the case of the modified-states solver. 
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Fig. 21. — Two-dimensional simulations of a perturbed nova precursor atmosphere, zoomed 
in to near the interface, evolved 0.20 s. Plotted is temperature with velocity vectors. A 
white contour shows the highest temperature region (T > 4.01 x 10 7 K), and a black contour 
shows the interface between the C/O white dwarf and the accreted material. Above are the 
simulation results without adjusting the profile to HSE in this code, with reflecting boundary 
conditions, and without the modified-states PPM. Below are results with the second-order 
differencing to bring the model to HSE, second-order constant-temperature HSE boundary 
conditions, and using the modified-states PPM. Note the difference in velocity magnitudes, 
spurious compression of the bottom layer, and resulting compressional heating, which wipes 
out the real physical effect under consideration, the initiation of convection. 



